Is the Yb 2 Ti 2 07 pyrochlore a quantum spin ice? 
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We use numerical linked cluster (NLC) expansions to compute the specific heat, C(T), and en- 
tropy, S(T), of a quantum spin ice model of Yb2Ti2 07 using anisotropic exchange interactions 
recently determined from inelastic neutron scattering measurements and find good agreement with 
experimental calorimetric data. In the perturbative weak quantum regime, this model has a ferri- 
magnetic ordered ground state, with two peaks in C(T): a Schottky anomaly signalling the para- 
magnetic to spin ice crossover followed at lower temperature by a sharp peak accompanying a first 
order phase transition to the ferrimagnetic state. We suggest that the two C(T) features observed 
in Yb2Ti2 07 are associated with the same physics. Spin excitations in this regime consist of weakly 
confined spinon-antispinon pairs. We suggest that conventional ground state with exotic quantum 
dynamics will prove a prevalent characteristic of many real quantum spin ice materials. 

PACS numbers: 74.70.-b,75.10.Jm,75.40.Gb,75.30.Ds 
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The experimental search for quantum spin liquids 
(QSLs), magnetic systems disordered by large quan- 
tum fluctuations, has remained unabated for over twenty 
years [if. One direction that is rapidly gathering mo- 
mentum is the search for QSLs among materials that are 
close relatives to spin ice systems [2[ , but with additional 
quantum fluctuations, or quantum spin ice |3j, |j]. 

Spin ices are found among insulating pyrochlore oxides, 
such as R 2 M 2 7 (R=Ho, Dy; M=Ti, Sn) @. In these 
compounds, the magnetic R rare earth ions sit on a lattice 
of corner-sharing tetrahedra, experiencing a large single- 
ion anisotropy forcing the magnetic moment to point 
strictly "in" or "out" of the two tetrahedra it joins (see. 
Fig. la). Consequently, the direction of a moment can be 
described by a classical Ising spin [2j. In these materials, 
the combination of nearest-neighbor exchange and long- 
range magnetostatic dipolar interactions lead to an expo- 
nentially large number of low-energy states characterized 
by two spins pointing in and two spins pointing out on 
each tetrahedron (see Fig. la). This energetic constraint 
is equivalent to the Bernal-Fowler ice rule which gives 
water ice a residual entropy «Sp ~ ^3(0) ln(3/2) per pro- 
ton, estimated by Pauling [6 1 and in good agreement with 
experiments on water ice [2|. Since they share the same 
"ice-rule", the (Ho,Dy)2(Ti,Sn)207 pyrochlores also pos- 
sess a residual low-temperature Pauling entropy Sp |8(, 
hence the name spin ice. The spin ice state is not thermo- 
dynamically distinct from the paramagnetic phase. Yet, 
because of the ice-rules, it is a strongly correlated state 
of matter - a classical spin liquid of sorts [l|, [2J . 

For infinite Ising anisotropy, quantum effects are ab- 
sent |2j. However, these can be restored when consid- 
ering the realistic situation of finite anisotropy. In two 
closely related papers, Hermele et al. |9| and Castro-Neto 
et al. [10| considered effective spins one-half on a py- 




FIG. 1: (a) Two neighboring tetrahedra with spins in their 
two-in/two-out ground state, (b) spinon/antispinon pair, (c) 
spinon/antispinon pair separated by a (green) string of mis- 
aligned spins in the pyrochlore lattice. 



rochlore lattice where the highly degenerate classical spin 
ice state is promoted via quantum fluctuations to a QSL 
with fascinating properties. This QSL is described by 
a compact lattice quantum electrodynamics (QED) -like 
theory. In this QSL state inherited from the parent clas- 
sical spin ice, the ice-rules amount to a divergence-free 
coarse-grained fictitious electric field whose sources are 
deconfincd spinons while the sources of the canonically 
conjugate field are deconfincd monopoles |ll|, along with 
a gauge boson ("artificial photon"). 

Recent numerical studies have found evidence that 



QED-like phenomena may be at play in some minimal 
quantum spin ice (QSI) lattice models [13j - but does 
the QSI picture apply to real materials? Also, should 
a QSI state be solely defined by whether or not a QSL 
state is realized? While a QSI picture has been sug- 
gested relevant to the QSL behavior in Tb2Ti207 [3J 
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intense experimental |14421| | and the- 
2ll426j interest has recently turned to 
Yb 2 Ti 2 07 (YbTO), which has been argued to be on the 
verge of realizing a QSL originating from QSI physics. In 
fact, the combination of (i) an unexplained transition at 



T c ~ 0.24 K [1J, [27(, (ii) the controversial evidence for 
long-range order below T c [28j,|29[ and (iii) the high sensi- 
tivity of the low-temperature (T < 300 mK) behavior to 
sample preparation conditions [13, 123] are all tantalizing 
evidence that YbTO has a fragile and perhaps uncon- 
ventional ground state. Thus, explaining YbTO is a key 
milestone in the study of QSI in a materials context. 

In this paper, we first use the numerical linked clus- 
ter (NLC) method [33, |31| to calculate the heat capac- 
ity, C(T), and entropy, S(T), of a microscopic model 
for YbTO with exchange parameters, {J e }, taken from 
Ref. [l8(. This calculation, which converges down to 
about 1 K, agrees well with experiments. It demonstrates 
that YbTO is indeed a spin-half, anisotropic exchange 
model, with {J e } determined from magnon energies in 
the strong-field polarized paramagnet regime 18 1. Our 
work suggests that a two-peaked C(T) structure is natu- 
ral in YbTO and should be present in the best ( "quality" ) 
samples [l3, |23J. Below the higher temperature C(T) 
hump near 2 K, the system has a residual S(T) com- 
parable to Sp, but without a clean S(T) w Sp plateau 
developing upon cooling. We propose that the lower tem- 
perature sharp peak in C(T) is associated with a strongly 
first order transition to a ferrimagnetic state. Such a 
behavior is indeed found in our study when the quan- 
tum (non-Ising) exchanges are small. Finally, we argue 
that despite a conventional ground state, the spin excita- 
tions consist of spinon/antispinon pairs connected with 
(Dirac-like 12|) strings of reversed spins, whose confine- 
ment length l s diverges in the limit of small quantum 
exchanges. We propose that these excitations should ul- 
timately form the basis for describing what we expect to 
be highly unconventional inelastic neutron spectra [261 ] . 

Model & Method - The anisotropic exchange QSI 
model is defined by the nearest-neighbor Hamiltonian 



18L l25l ] on the pyrochlore lattice 

H Q si = Y, i J ** S i S j - X J±(S+SJ + S^S+) 
+\J±±[~fijS?Sf + lijSrSr] 

+xj z± [(S!(q j s+ + QST) + 1 o j}}. (1) 

is a 4 x 4 complex unimodular matrix, and C = — 7* 
[181 ] . The z quantization axis is along the local [111] direc- 
tion, and ± refers to the two orthogonal local directions. 
We take A = 1, except when stated otherwise. 

Recently Ross et al. [18j used inelastic neutron scat- 
tering data in high fields to deduce the { J e } exchange pa- 
rameters for YbTO: J zz = 0.166±0.04, J± = 0.05±0.01, 
J ±± = 0.05 ± 0.01, and J z± = -0.14 ± 0.01, all in mcV. 
These parameters have also been determined through an 



analysis of the zero-field energy-integrated paramagnetic 
neutron scattering [l7j,l21|, but the values of the {J e } pa- 
rameters disagree significantly - an issue that we address 
in the supplementary material [32j . 

NLC expansions provide a controlled way of calculat- 
ing macroscopic properties of a thermodynamic system 
33,[3l|. By summing up contributions from clusters upto 
some size, one can obtain properties in the thermody- 
namic limit, which include all terms in high temperature 
expansions upto some order. Furthermore, since the con- 
tributions of the clusters are entirely included for all tem- 
peratures, all short distance physics is fully incorporated, 
and thus can converge down to lower temperatures than a 
(high-temperature, T) series expansion [30( in 1/T. NLC 
is particularly suited to the study of spin ice systems. It 
was recently shown that for classical spin ice models, just 
first order NLC based on a single tetrahedron, gives C{T) 
and S(T) for all T within a few percent accuracy [33j. 

Here, we calculate the thermodynamic properties of 
the exchange QSI model of Eq. ([6]) using tetrahedra- 
based NLC upto 4 th order [321 . Euler extrapolations J34J 
are used to eliminate some alternating pieces in the ex- 
pansion, which further improves the convergence of the 
calculations to lower T. In zero field, there is only one 
cluster in each of the first three orders, and three clusters 
in the fourth order [32|. The different g-tensor elements 
on different sites (expressed in a global frame) [2J] mean 
that many more clusters are needed for calculating field- 
dependent C(T), magnetization and susceptibility, and 
these will be presented elsewhere. 

Figure 2 shows C(T) calculated with different NLC 
orders. By 4 th order, there is good convergence to tem- 
peratures below the C(T) peak at ~ 2 K. Applying Eu- 
ler transformations |34| improves the convergence down 
to slightly below 1 K. The experimental data from Refs. 
271 ] . shown for comparison, agree well with the NLC re- 
sults. Here, we used the mean values of the {J e } from 
Ref. [l8[ and did not adjust any parameters. Given the 
variability in the experimental C(T) data from one group 
to another [l9ti2lll32l |. it does not seem useful at this time 
to search for {J c } parameters giving a better fit. This 
agreement shows that the { J e } parameters are not sub- 
stantially renormalized compared to the h igh (5 Tcsla) 
field values [1S|] . Using the {J c } of Refs. [17J, [2l( gives 
substantially different C(T) results [321 ] . 

Figure 3 shows S(T) calculated by NLC, together 
with the entropy obtained by integrating C (T)/ T data 
of Ref. [i3. We found the data from Ref. [27| ideally 
suited to perform this comparison [321 ] . The entropy con- 
verges to lower temperature slightly better than C(T) 
where, with Euler transformations, S(T) converges down 
to about 0.7 K, matching well with the experimental en- 
tropy values over the overlapping temperature range. 

Perturbative considerations - In order to better un- 
derstand the properties of this system, we turn to the 
perturbative regime A -C 1 in Eq. 1 18j,|25|. To second 
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FIG. 2: Specific heat, C(T), per mole of Yb for the model 
parameters in Ref. [18|]> m units of the Boltzmann constant 
&b, calculated via NLC (up to 4 th order NLC together with 
Euler extrapolations) are compared with experimental data 
for Yb 2 Ti20 7 . The black circles are data from Ref. [271 . 



order in A, only J z ±, by far the largest quantum term 
for YbTO, leads to a degeneracy- lifting classical poten- 
tial for different spin-ice configurations. It amounts to 
a fluctuation-induced ferromagnetic exchange constant 
,h = —3\ 2 J%±/J ZZ [25| between shortest distance spins 
on the same tetrahedral sublattice that share a neigh- 
bor [35|. It leads to the selection of a q = long-range 
ordered ground state in which all tetrahedra are in the 
same configuration and the spins develop a small ferro- 
magnetic moment along one of the (100) cubic directions. 
This q = ferrimagnet (FM) lacks the Coulombic physics 
originally present in the J zz -only spin ice model [36[. 

To calculate C{T) and S(T) in the perturbative regime 
at low T, we turn to classical loop Monte Carlo simula- 
tions [37| of the J3 



J zz model [32j. These reveal a 



very sharp lower temperature peak signalling a first or 
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der phase transition to a q = state (see Fig. S5 

Excited states in the perturbative regime: spinous and 
strings - A surprise of the perturbative treatment is that, 
while the ground state is classical, the spin-flip excita- 
tions remain non-trivial and of quantum nature. This is 
because, once a spin is flipped in a spin-ice state, creating 
a spinon/antispinon pair [11], the pair can hop through 
J z ± acting through first order degenerate perturbation 
theory. Thus, the dispersion in the excited state man- 
ifold is XJ Z ±, much larger than the dispersion within 
the low-energy manifold of spin ice states, which is only 

^ J z ±/Jz Z - 

A sketch of a spinon/antispinon pair is shown in Fig. 
lb and lc. Note that only spins inside the tetrahedron 
"already" containing spinons arc flippablc in first order 



FIG. 3: Entropy, S(T), per mole of Yb, in units &b following 
the methods described in the caption of Fig. 2. The black 
circles are obtained by integrating the data from Ref. [27[ 
excluding the nuclear (hyperfine) contribution. The Pauling 
entropy Sp ~ t? In | is shown as a horizontal line. The inset 
shows S(T) in the perturbative regime with Ja/J zz = —0.001. 
A clear plateau at S(T) » Sp is seen, followed at lower T by 
a precipitous drop of S(T) (i.e. latent heat) accompanying 
the transition to long range FM order |32l |. 



degenerate perturbation theory. Hence, the connecting 
string of misaligned spins can only fluctuate by higher 
order processes involving closed loops with alternating 
in-out spins [261 ] . Thus the renormalizcd string tension 
per unit length remains finite and of order J3. One can 
estimate the typical string length as the length, l s , at 
which the cost of the string becomes comparable to the 
derealization energy of the spinon/antispinon pair. The 
string energy per unit length goes as ~ J3 ~ A 2 , whereas 
the derealization energy (spinon bandwidth) goes as A. 
This leads to l s scaling as 1/A, which diverges as A — >• 0. 

A detailed theory of neutron scattering in this ferri- 
magnetic phase is not attempted here, but we anticipate 
it to follow the proposal of Ref. |26|. At temperatures 
above the transition to the q = long-range ordered 
state, the system explores the classical two-in/two-out 
spin ice states and should display singularities (pinch 
points, PPs) in neutron scattering [36[ rounded off by 
the finite density of thermally excited spinon/antispinon 
defects [IJ, [36| . While the system has thermally smeared 
PPs above the ferrimagnetic transition and no static PPs 
well below the transition, it may display some remnant 
of PPs in the spin dynamics at higher energies. These 
interesting issues deserve further attention. 

Beyond the A <C 1 regime - Why is the transition 
temperature of YbTO so low? As discussed by Ross et 
al. [l8j . the low T peak in C{T) is at a temperature lower 
than mean-field theory by an order of magnitude. Com- 




FIG. 4: Monopole defect density, p(T), calculated using NLC, 
shown down to a temperature where 3rd and 4th order Euler 
Transforms agree. Here, quantum exchanges are scaled with 
respect to YbTO parameters by different values of A. 



paring C(T) for the quantum model with different A with 
the corresponding classical model with the perturbative 
J3/ ' Jzz value provides a hint of the reason why [321 ] - It 
shows that, in the classical model, the long-range order 
keeps steadily moving up with increased J3, even beyond 
the short-range order C(T) peak. In contrast, the quan- 
tum systems, with different A continue to display a short- 
range order C(T) peak and presumably long-range order 
only occurs at a much lower T. Perturbative considera- 
tions here have an analogy with strong coupling studies 
of Mott physics in the Hubbard model, where the Neel 
temperature first increases with t as t 2 /U but then be- 
gins decreasing when the system moves away from the 
perturbative small t/U regime. We propose that a sim- 
ilar non-monotonic T c arises in this QSI model due to 
enhanced quantum fluctuations. 

Another argument for a reduced T c comes from 
considering the temperature dependence of the defect 
(spinon/antispinon) monopole density, p(T), as calcu- 
lated by NLC (see Fig. 4 and Figs. S3 and S4 [H). To 
illustrate the point, we show the behavior for several dif- 
ferent A values. Convergence increases to lower T, with 
decreasing A, as expected. One finds that as T drops 
below the hump in C(T), p(T) displays a plateau-like 
region, whose value increases steadily with increasing A. 
This indicates that the states within the spin-ice man- 
ifold develop large spinon/antispinon spectral weight, 
thus strongly renormalizing all low energy scales and, 
presumably, leading to reduced T c . 

Discussion: What constitutes an exchange QSI? - 
We suggest that a double-peaked C(T) with an en- 
tropy between the peaks comparable to 5p is the hall- 
mark of an exchange quantum spin ice (QSI). How- 
ever, one is unlikely to find an exact plateau at S(T) « 



Sp outside the perturbative (small A) regime. Such a 
double-peaked structure and quasi-separation of the en- 
ergy/temperature scales associated with short and long- 
range physics has also been suggested for other systems 
where quantum spin liquid physics may apply |38| . 

According to the gauge mean-field theory of Ref. [25| , 
at low temperature below which short-range spin ice cor- 
relations develop, a system may exhibit either a conven- 
tional ferrimagnetic (FM) order, a Coulombic ferromag- 
net (CFM) or a full-blown quantum spin-liquid (QSL), 
depending on its quantum exchange parameters. The 
largest quantum exchange terms in YbTO is J z ±, which 
favors the FM state, which we believe is the origin of 
the 0.24 K transition in the best samples [28[. It re- 
mains to be seen if there are real materials for which J± , 
which favors the QSL [3, llfj, [25[ , is the dominant quan- 
tum term. Nevertheless, even when the ground state is 
FM, the excitations remain highly exotic, consisting of 
spinon-antispinon pairs separated by long strings. This 
non-trivial feature is derived from the underlying spin-ice 
physics. Finally, as one notes that J z ± is strictly zero for 
non-Kramers ions (e.g. Pr, Tb) and that virtual crystal 
field excitations [3J in Tb-based pyrochlores are a fun- 
damentally different pathway from anisotropic superex- 
change [4J to generate anisotropic { J e } couplings between 
effective spins one-half [3j, |4[ , the prospect to ultimatel 
find a QSI-bascd QSL among rare-earth pyrochlores 
is perhaps promising. 
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DMR-1004231, the NSERC of Canada and the Canada 
Research Chair program (M.G., Tier 1). We acknowledge 
very useful discussions with B. Javanparast, K. Ross and 
J. Thompson. Wc thank P. D almas de Reotier for pro- 
viding specific heat data of Ref. [L9( . 



Supplementary Material 

This supplement provides the reader with further ma- 
terial to assist with some of the technical materials of the 
main part paper 



Numerical Linked Cluster Method 

For the proposed QSI Hamiltonian 1181 ■ the numeri- 
cal linked cluster (NLC) method [30j, |3l| gives reliable 
quantitative properties of the system in the thermody- 
namic limit down to some temperature by developing an 
expansion in connected tetrahedra that embed in the py- 
rochlore lattice. For each cluster, we perform an exact 
diagonalization (ED) and calculate physical quantities 
from the resulting spectrum and states. Once a prop- 
erty is calculated, the properties of all subclusters arc 
subtracted to get the weight of the cluster c denoted as 



W(c). In the thermodynamic limit, an extensive prop- 
erty, P is expressed as 



P/N = J^L(c)xW(c), 



(2) 



where L c is the count of the cluster, per lattice site. 

We consider all clusters up to four tetrahedra, the 
largest diagonalization being a 13-site system. All states 
are required to calculate the partition function and ther- 
modynamic quantities presented below. The particular 
clusters to fourth order in our expansion are shown in 
Figure SI. 



Computational Requirements 

NLC using the tetrahedral basis requires exact diago- 
nalization of increasingly large tetrahedral clusters. Us- 
ing modern hardware and freely-available linear algebra 
routines, diagonalizations for clusters of one tetrahedron 
(four sites) and two tetrahedra (seven sites) could be done 
in less than a second, while the three-tetrahedron (10- 
site) cluster still required less than 10 seconds. Comput- 
ing only the spectrum for a single four-tetrahedron (13- 
site) cluster required about 1200 seconds and more than 1 
GB of memory, while generating the full set of eigenstates 
required approximately 8 GB of memory. Note that the 
Hamiltonian of an N-site cluster is a 2 N x 2 N complex 
Hermitian matrix. Exact diagonalizations of larger sys- 
tems are, in practice, limited by memory requirements. 
The next order calculation will have 3 more sites and the 
memory requirement will grow by a factor of 64. 




FIG. 5: SI: Clusters used for the zero-field NLC expansion 
in the tetrahedral basis, up to fourth order. Each graph is 
accompanied by its lattice constant L. 



series. In our case, where direct terms are available 
to fourth order, we begin the Euler transform after the 
second order, so that the third and fourth order Eulcr- 
transformed property estimates are 

°3,E = Sq + Si + S-2 + ySs, 



R 



4.E 



P 



S3 + S4 



3,E 



(5) 



Euler Summation 



NLC generates a sequence of property estimates {P n } 
with increasing order n, where P n = J27=i &i an d Si is 
some physical quantity calculated at the ith. order. When 
such a sequence is found to alternate, its convergence can 
be improved by Euler Transformation |34| . In general, 
given alternating terms Si = (— l) l Ui, the Euler Trans- 
form method amounts to estimates, 



Uq — U\ + u-z — . . . — u 



«-i+E^- [ASu - ] 



(3) 



where A is the forward difference operator 
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A 2 u„ 
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U n +2 


- 2u n+ i +u n , 


A 3 u„ 
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U n +3 


- 3u n+ 2 + 3w„+i 



(4) 



Usually, a small number of terms are computed directly, 
and the Euler transformation is applied to rest of the 



Various Hamiltonians and perturbative limit 



We use the notation of Ross et al. 
quantum spin ice Hamiltonian as 



1181 and define the 



%QSI 



<i,j> 



{J zz S*S*-J±(StSr + S7Sf) 



+j±±hi j stsf + 7 * j s-s-] 

+J*±[(S? (Ca S^ + Ctj S~) + i^j}}. (6) 

The parameters for Yb2Ti20y determined by fitting from 
high-field inelastic neutron (magnon) spectra in Ref. [18[ 
are, measured in meV, J zz = 0.166 ± 0.04, J± = 0.05 ± 
0.01, J±± = 0.05 ± 0.01, and J z ± = -0.14 ± 0.01. Two 
other sets of parameter estimates for Yb2Ti207 were de- 
termined by fitting the diffused (energy-integrated) neu- 
tron scattering using the random phase approximation 
(RP A) [fR EH. The values obtained by Thompson et 
al. [13 are: J zz = 0.023, J± = 0.038, J ±± = 0.007, 
and J z ± = —0.040, while those obtained by Chang et 
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FIG. 6: S2: Molar heat capacity for YbTO reported by Blote 
et al. [27J and by Yaouanc et al. [19l | compared with cal- 
culated values using exchange parameters from Ross et al. 
(Ross-E3,E4) [131, Thompson et al. (Thompson-E3,E4) [l7[ 
and Chang et al. (Chang-E3,E4) [5l[. Third (E3) and Fourth 
(E4) order Euler Transforms of the NLC results using the pa- 
rameters are shown. 



al. [2l| are J ^ = °- 059 ' J ± = O- 023 ^ J ±± = °- 006 > 
and J z ± = —0.029. In all cases, the values of the { J e } 
exchange parameters are given in meV. The calculated 
heat capacity for all these parameters, together with the 
experimental data on Yb2Ti20y from difference groups 
19l |20(, are shown in Fig. S2. It is clear that the latter 
two parametrizations by Thompson et al. and Chang et 
al. do not give a good description of the heat capacity of 
the material. It is not clear at this time why RPA calcu- 
lations find such { J e ) parameters compared to high- field 
paramagnon spectra [20j . This problem warrants further 
attention. 

In order to explore to what extent quantum mechanical 
effects are at play in Hqsi, we introduce a Hamiltonian 
with rescaled quantum terms as 



H\ — Ho + XJii, 



(7) 



where Ho is the classical spin-ice Hamiltonian consisting 
of J zz terms only, while all other terms are included in 
Hi. The value A = 1 corresponds to the parameters of 
Ross et q?.|18j In the perturbative regime (A < 1), this 
model maps on to a J\ — J3 model with J\ = J zz and 
J3 = — 3A J z ±/J zz . 

Specific heat and entropy of the system with different 
values of A in 4th order Euler Transform, down to a tem- 
perature where 3rd and 4th order Euler Transforms agree 
with each other are shown in Fig. S3 and Fig. S4. Heat 
capacity of the perturbative classical J\ — J3 model, cal- 
culatcd by classical loop Monte Carlo simulations [37[ is 



FIG. 7: S3: Heat capacity where quantum terms are scaled 
by A. 
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FIG. 8: S4: Entropy of the system, when quantum terms are 
scaled by A. The orange line is the Pauling entropy S v . 



shown in Fig. S5. Note that while the models with differ- 
ent A always have a short-range order peak, in the ,J\ — J3 
model, long-range order temperature increases well past 
the short-range order peak with increasing J3/J1. 



Comparison of the experimental entropy vs NLC 
results 



The entropy difference, SiT-i) — S{T\) between two 
temperatures T\ and Ti can be obtained by integrating 



C(T)/T between those two temperatures: 



Monte Carlo Simulation of the J zz — Js Model 



5(T 2 ) - 5(Ti) = f 2< ^P-dT 



ii 



T 



The number of experimental specific heat, C(T), re- 
sults on Yb gTig O? has rapidly accumulated over the past 
year or so 19|-|21[. Most of these data are somewhat 
problematic in wanting to assess whether those thermo- 
dynamic data hide spin ice phenomenology, associated 
with a rapid diminution of spinon/antispinon excitation 
and the concurrent C(T) hump at a temperature ~2 K 
as we now discuss. 

All of the published C(T) data [H-ED, Ez| do not go 
to sufficiently high temperature to extract reliably the 
limiting C(T) oc 1/T 2 high temperature behaviour that 
would allow one to determine the residual magnetic en- 
tropy by integrating C(T)/T upon decreasing T starting 
from the infinite fee hi(2) value. One must therefore in- 
tegrate C(T)/T from low temperature, and assume an 
entropy value, iSi ow at some reference (low) temperature, 
Xi ow . The apparent large amount of residual entropy be- 
low ~ 0.2 K in the single crystal samples of Refs. [19l — 
2l| make difficult ascribing a reasonable value to S\ ovl . 
This problem is further compounded by the rising low- 
temperature nuclear contribution to the total specific 
heat below about 0.1 K. The very sharp 1st order transi- 
tion seen in powder powder sample of Ref . |20( , without 
a precise measurement of the associated latent heat also 
make difficult using those data for comparison of experi- 
mental entropy with the S(T) calculated by NLC. On the 
otherhand, the data of Blote et al. [27J seem the most ade- 
quate for comparison with NLC: there is a sharp specific 
heat peak at T c ~ 0.24 K with sufficient temperature 
resolution that allows integration of C{T)/T over the 
peak without concern about an associated latent heat. 
The C{T) data are dropping rapidly below T c , suggest- 
ing the opening of an excitation gap, ultimately reaching 
a low- value that is limited by the "high temperature tail" 
(T ~ 0.1 K) of the nuclear contribution. Using the data 
from Ref. [27| , we thus assume that the magnetic part of 
the specific heat is zero at T = 0.1 K, and integrate up- 
ward (increasing temperature) C(T)/T up to the highest 
temperature point available from those data (~ 3.5 K). 
This results in the data (filled black circles in Fig. 3 in 
the body of the paper). 

It would be highly desirable to repeat this procedure 
from the C(T) data of Refs. 19h21| which show a sharp 
peak, but including (magnetic specific heat) data for T 
up to 20 K where the limiting high-temperature regime 
C(T) w ^ + ^ can be fitted and compared with NLC, 
along with measurements of the magnetic entropy, S(T). 

limit available dataThe data from Working from the 
reasonable presumption high temperature 



In the perturbative regime of the QSI, we consider the 
effective Hamiltonian 



H = 



E 



■/, 



E 



■ho id i 



(8) 



where a = ±1 are the Ising variables. (. . .) denotes 
the sum over the nearest neighbors, (. . .)' denotes the 
sum over the third nearest neighbors which share a near- 
est neighbour. Distance-wise there exists another type 
of third nearest neighbors which do not share a nearest 
neighbor. For any given spin, there are six third near- 
est neighbors for both types. Antifcrromagnctic J zz > 
drives the spin ice formation in the classical spin ice sys- 
tem, and a small fluctuation-induced ferromagnetic ex- 
change Jj, = —3J Z± /J ZZ < favors the q = ordering 
within the spin ice manifold, i.e., all tctrahedra on the 
same primitive FCC lattice have the same one of the six 
spin ice states. 

Monte Carlo simulations are performed using the 
Metropolis algorithm. Single spin flip updates are used 
along with the non-local loop algorithm [37j . which re- 
stores the ergodicity of the system once it is frozen into 
the spin ice states. Systems of 128 spins are simulated 
in a cubic box with periodic boundary conditions. Up 
to about 78,000 Monte Carlo steps per spin are used in 
equilibrating the system at a given temperature, with the 
same number of steps in data sampling. To investigate 
the calorimetric quantities, fluctuations of the energy are 
recorded to give the heat capacity: 



C = 



< E 2 > -< E >' 
k^f 2 



(9) 



Calculation of Monopole density 

The defect (spinon/antispinon) monopole number 
M(T), for a cluster, is evaluated as 

M(T) = tr(fhe- /3fl )/Z 



We-f'-Mktfm,, 



(10) 



where mu is the monopole count in the local S z basis 
state \k). This count is a sum over all the tetrahedra in 
a cluster, mk = X)j m ki, where 

2 all in/out, 
mki = { 1 three in/out and one out/in, (11) 

two in and two out. 
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FIG. 9: S5: Heat capacity of the classical J zz — Jz model, 
with different J zz /Ji ratios. 



The monopole density p(T) is defined as number of 
monopolcs present per site, giving 



p(T) = M(T)/N. 



(12) 
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